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Abstract 

Dryland ecosystems commonly exhibit periodic bands of vegetation, thought 
to form due to competition between individual plants for heterogeneously 
distributed water. In this paper, we develop a Fourier method for locally 
identifying the pattern wavenumber and orientation, and apply it to aerial 
images from a region of vegetation patterning near Fort Stockton, Texas. We 
find that the local pattern wavelength and orientation are typically coherent, 
but exhibit both rapid and gradual variation driven by changes in hillslope 
gradient and orientation, the potential for water accumulation, or soil type. 
Endogenous pattern dynamics, when simulated for spatially homogeneous 
topographic and vegetation conditions, predict pattern properties that are much 
less variable than the orientation and wavelength observed in natural systems. 
Our local pattern analysis, combined with ancillary datasets describing soil 
and topographic variation, highlights a largely unexplored correlation between 
soil depth, pattern coherence, vegetation cover and pattern wavelength. It also, 
surprisingly, suggests that downslope accumulation of water may play a role in 
changing vegetation pattern properties. 



1. Introduction 



Pattern formation occurs in numerous ecological and biological systems, where it 



has been linked to reaction-diffusion (Turing) type instabi lities dRietkerk et al. 
|2002| ), hydrodynamic instabilities ( [Thompson & Daniels[ |2010| ), and potentia- 



t e variational) dynamics that maximize or at le ast increase ecological productivity 
( [Lefever et al.\ |2009[ |Pringle et al\ |2010| ). Patterns have be en observed in 
drylan d vegetation (|B orgogn o et <3/.[|2009| ), in bogs and wetlands ([Rietkerk et al. \ 
2004Z?t jLarsen & Harvey [ |201 in mussel beds ( [Van de Koppel et al \ |2008L 



termite mounds ([Pringle fl/.[|2010| ) and other systems ( |Rietkerk et al.\\20QAcfy . 
In many cases ecological patterns form an intermediate realization between 
environmental states in which the entire landscape is either colonized or bare. As 
such they often indicate the presence of bistable states, which are characterized 



by the potential for critical and locally irreversible transitions (Schef fer et al. 
2009 ). Observations and theoretical work both suggest that transitions between 



vegetated and desertified (unvegetated) conditions in patterned systems are 
preceded by striking changes in the morphology of the vegetation patterns, which 



thus a ct as an early-warning sign of deteriorating ecosystem health ( |Kefi et al. 
2007). Understanding the controls on the morphology of vegetation patterns 



therefore has practical interest in terms of ensuring that observed changes are 
interpreted correctly. 

Vegetation pattern formation in dryland ecosystems is a global pheno menon, 
ranging from random distributions of bare soil and vegetation canopies ( |Caylor| 
\et al.\ 2004| ) to highly organized spatial distributions with identifiable length- 
scales and orientations (HilleRisLa mbers et al.\ |2001[ |Rietkerk et al.\ |2002; 
Borgogno et al.\ |2009|). Int ermediate cases such as power-law (scale-free) 



cluste ring ( [Scholes & Archer] |1997^ |Scanlon et al.\ |2007[ |von H ardenb erg et al. 
2010) , and dendritic structures in which vegetation co ncentrates along drainage 



lines ( |McGrath et al.\ |2012| |Thompson et al.\ |2011a| ) are also observed. All of 
these morphologies can be related to the presence, strength and directionality of 
positive feedbacks that concentrate resources (such as nutrients, soil carbon anc 



water) that sustain plant life in a localized region near the plants (Ravi et al. 



2008 



2005 



Greenej [T992| |Galle et al\ [T999| |Harman et al\ |20T2l |Fuigdefa bregas 



Schlesing er et al.\ |1990| ). These feedbacks have lead to the moniker 
'ecosystem engineers' being applied to perennial plants in dryland ecosystems: 
they create the conditions necessary for their own survival ( |Gilad et aL\ |2004 



Yizhaq et al. 2005) 



The formation of periodic vegetation patterns is strongly linked to the 
coupling of (i) water redistribution to the vicinity of plants, an d (ii) competition 
between individual plants for access to this water res ource ([Borgogno et al. 



2009] |Bromley et ~aL\ [T997t |Lefever et al\ [2009t |Seghieri et al.\ |1997| ). 
The dynamics of patterned vegetation systems has formed a focus of 
ecohydrological and nonlinear dynamics research, motivated by the increasingly 
well-demonstrated connection between pattern morphology and desertification 
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et 



sk flRietkerk et al\ [20021 |Kefi et al\ [20071 jBarbier et al\ [20081 jDeblauwe 



\et al\ 2011); and by the inherently interesting n onlinear dynamics of these 
systems fLefever et al\ |2009| |Meron et al\ |2004| ). Field and theoretical work 
has identified the important roles of vegetatio n dynamics (|Kefi et al\ |2007| ), 
seed dispersa l processes ([ Thomps on et al. [2008[ [Thompson & Katul||2009| ), root 
morphology (|Barbier et fl/.[|2008[|Lefever et fl/.[|2009|), surface flow dynam ics 
( |Thompson et al\ |201 and climate feedbacks ( |Konings et al\ |2011[ ) in 
modulating pattern morphology. 

While our understanding of vegetation pattern dynamics is improving, there 
remain undeniable differences between simulated vegetation patterns and the 
natural ones observed in the field. Natural patterns are characterized by a 
degree of disorder and heterogeneity on multiple scales that is not reproduced 
in idealized models. Disorder can arise intrinsically through the existence of a 
range of stable wavelengths and orientations, with transitions in space and time 
occurring through local pattern pattern defects ( [Cross & Hohenbergl |1993| ). In 
addition, the pattern wavelength and orie ntation can be stro ngly influenced by 
either wavelength-scale in homogeneities ( |Lowe et al\ |1983| ) or larger gradients 
( |Eckhaus & Kuske[|1997|) i n the und erlying system, as well as the presence of 
boundaries ( |Hu et a/.[|1993| ) or noise ( |Lindner[|2004| ). 

In natural systems, this means that disorder observed in vegetation patterns 
could either reflect intrinsic features of the pattern forming process {endogenous 
effects), or could reflect spatial c hanges in soil structure and local topography 
{exogenous effects). For instance, Thom pson et <z/.| ( [2008| ) explored whether the 
unrealistically smooth nature of many models of vegetation biomass distribution 
was an artifact of representing seed dispersal as a diffusive process. Greater 
fidelity between modeled and observed (disordered) vegetation patterns was 
achieved by representing plant population migration with a seed dispersal kernel. 
More simply, theoretical treatments of vegetation patterning usually neglect 
variations in soils and topography, or impose periodic boundary-conditions that 
remove the differences in water availability between the top and bottom of a 
hillslope. Such simplifications will inevitably lead to idealized representations 
of the vegetation response and obscure interactions between the intrinsic pattern 
structure and the spatial structure of the landscape. For example, McGrath 
\et al\ ( [201 2| ) demonstrated that the orientation and wavelength of vegetation 
patterns modeled with realistic boundary conditions changed between the top and 
bottom of the hillslope. In control simulations with periodic boundary conditions, 
vegetation bands formed with a single wavelength and were orientated at 90° to 
the direction of water flow, in agreement with other modeling studies. When 
realistic boundary conditions were imposed, however, the bands near the top of 
the hillslope curved to lie perpendicular to the ridgeline. Similarly, the effects 
of spatial heterogeneity in soil properties on pattern formation have not been 
widely explored (although see [Thompson et al\ ( [20 08)). A number of theoretical 
studies indicate that local biomass, band properties and band coherence should 
vary with changing minimum and maximum infiltration capacities and soil 
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properties ( |Ursino & Contarinij |2006[ |Thiery et al\ |1995| ), and there are some 
tantalizing hints that subsurface features, such as the ironstone underlying tiger 
bush in Niger, calcrete hardpans underlying banded patterns in Texas, and 
silcrete hardpans underlying mulga bands in Australia ( |McDonald et al\ |2008 ; 
White} [T970t|Mabbutt & Fanning[|1987| ), may have an association with patterned 



vegetation. 

Predictions about the interaction of vegetation patterns with changing soil 
or topography can be investigated using remotely sensed datasets, an approach 
with a long and growing history. Vegetation patterns in Africa were first 
observed from light aircraft flights (Macfadyen} |1950[ |Worrall[ |1960| ). Initial 



analyses of the patterns demonstrated their spatial regularity on the basi s of the 
two-dim ensional Fourier powe r spectrum ( [Couteron & Lejeunel |2001| ). More 



recently, Deblau we et a/.| ( |20TT] ) undertook large scale analyses of morphological 
trends in vegetation patterns in the Sudan by linking several remotely sensed 
datasets: surface imagery from the System for Earth Observation (SPOT) 
satellite, topographic data from the Shuttle Radar Topography Mission (STRM), 
and rainfall data from the Tropical Rainfall Measuring Mission (TRMM). This 
allowed analysis of pattern morphology at 400 m resolution over multiple square 
kilometers. Although higher resolution datasets are available, they have g enerally 



ither been analyzed only over relatively small spatial scales 1 km 2 ) ([Barbier 



\et al\ [2008 ) or in vestigated from the pe rspective of identifying morphological 



change over time ( [Deblauwe et a/.[|?012| ). 

A key open question is therefore to characterize small-scale irregularities 
in vegetation patterns, and to classify them based on whether they arise due to 
either endogenous dynamics or variation in exogenous features imposed by the 
landscape. Understanding the implications of such variation on the resilience 
and stability of the ecosystems has important consequences for identifying and 
preventing potential desertification. 

In this paper, we report on the development of methods suitable for 
quantifying local patterns within high-resolution aerial photography, and for 
relating those features to ancillary datasets describing soil and topographic 



variation. Our site, located near Fort Stockton, Teaxs (see McDonald et al. 
( |2008| )) was selected because it exhibits vegetation bands and also offers several 
characteristics that facilitate studies of covariation. First, the aerial photography 
covers a large region at a resolution (0.5 and 1 m) comparable to the diameter 
of perennial vegetation canopies, allowing for the observation of individual 
plants. Second, high resolution (10 m spatial, 15 cm vertical) elevation datasets 
are available from the US National Elevation Dataset, allowing changes in the 
orientation and gradient of the hillslope to be mapped on scales that are much 
less than the pattern wavelength. Finally, the study area has been mapped as 
part of the SSURGO national soils database, and contains considerable variation 
in soil type. The availability of these datasets provides a valuable opportunity 
to investigate correlations between vegetation patterning and soil characteristics 
over tens of square kilometers. 
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Our analysis depends on these spatial datasets, and is subject to their inherent 
limitations. These limitations include the resolution, spatial artifacts, and the risk 
of spurious correlation, given that vegetation features are often considered when 

135 mapping soil boundaries. To minimize the effects of issues, we focus on two 
large-scale hypotheses: (1) The wavelength and orientation of the vegetation 
pattern are locally coherent but exhibit both rapid and gradual variation; and (2) 
The variability in vegetation pattern features will correlate with soil and elevation 
features in predictable ways. We note that by focusing on soil and topographic 

140 features, we inherently assume that spatial variations in pattern morphology 
arise due to spatial heterogeneity in landscape properties, rather than due to the 
nonlinear dynamics of the pattern forming process itself. In patterns far from 
threshold, features such as defects, dislocations, grain boundaries, and boundary 
conditions can result in spatial heterogeneity in pattern properties, even when 

145 all other fields are homogeneous ( [Cross & Greensidel|2009| ). To control for this 
possibility, we also briefly address the following null hypothesis: (3) Large-scale 
variations in pattern properties can be explained by the nonlinear interactions 
associated with vegetation pattern formation. 

We test the major hypotheses through the development of a localized 

150 Fourier technique to identify pattern wavenumber and direction. To address 
Hypothesis 1, this technique was applied to high (0.5 m) resolution aerial 
photography. To address Hypothesis 2, the technique was applied to a 188 
km 2 area, allowing local pattern properties to be correlated to the site soil and 
topographic properties. Hypothesis 3 was tested by identifying a subset of the 

155 0.5 m resolution image containing banded patterns that were close to ideal (i.e. 
reproducible by a model). We applied the Fourier analysis technique to both 
observed and simulated patterns, and compared variability in the wavenumber 
and direction fields in the modeled and observed patterns. 

We find that our two major hypotheses are satisfied. Local patterns are 

leo oriented in the same direction as the topographic slope and the pattern 
wavelength decreases for for steeper gradients. Deviations from these trends are 
associated with the presence of ridges, stream channels, anthropogenic features 
or changes in soil type. Different soil types within the study area determine the 
pattern boundaries and the pattern morphology: shallow soils are associated with 

165 highly coherent, shorter wavelength patterns, and deep soils with patterns that 
become incoherent and increase in wavelength near stream channels. Finally, 
the modeling test validated our decision to focus on landscape and soil features. 
Even in the most uniform region of the observed patterns, the variability in the 
real pattern properties exceeded that which could be retained in the steady state 

170 solution of a physical model that closely reproduced the mean pattern properties. 
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2. Methods 



(a) Study Site and Data 

The study site is a 188 km 2 area located approxim ately 30 km NW of Fort 
Stockton, Tex as (coordinates: 31°05' N, 103°03' W) ( [National Climatic Data 
[Center| |2010| ). The climate is hot and dry, receiving 370 mm annual rainfall 
on average, mean summer maximum temperatures of approximately 35°C and 
winter minima near freezing. The site is part of a large cattle ranch and is used for 
grazing. Dominant vegetation species include tarbush (Flourensia cenua), bunch 
and sod grasses (Aristida purpurea, Bouteloua curtipendula and Scleropogon 
brevifolius), and mixed mesquite (Prosops is glandulosa) and juniper (Juniperus 
pinchotti) brush ( |McDonald et al \ |2008| ). The site contains a striking spatial 



pattern consisting of bands of continuous vegetation cover lying over bare soil 
( [McDonald g7^|2008] ); see Fig. [I] 



High resolution imagery (0.5 m and 1 m pixels) were obtained from 
igital Globe, and from the National Agricultural Imaging Project (National 



Agri cultural Imaging Pro gram, 2 010[ ). We used the highest resolution images 
from Digital Globe for fine-scale analysis and a comparison between modeled 
and observed pattern morphology. We analyzed the NAIP images, which cover 
the whole area, to relate local pattern properties to soil and topographic features. 

We classified the image pixels as 'vegetation' or 'no vegetatio n' using a 
supervised classification based on total brightness (Rich ards} |1999| ). We used 



brightness because the perennial vegetation was not actively photosynthesizing 
when the photographs were taken, meaning that standard vegetation indices 
could not discriminate the vegetated locations. An example of the resulting 
binary image is shown in Fig. [T] Darker colors represent vegetation cover and 
the lighter colors bare soil. The insets show an original and classified image over 
a 260 x 260 m 2 window. 



(b) Fourier Windowing Method 

To quantitatively test our hypotheses, we developed a quasi-local technique 
to measure pattern wavelength A and pattern orientation for the binary images. 
The technique provides information about the local wavevector k = k x x + k y y, 
similar to that provided by short-time Fourier transforms or wavelet based 
appr oaches used in timeseries analysis (All en & Rabiner[ |1977[ |Daubechies 
199 0t. Local wa vevector s are useful for classifying convection patterns 



Heutmaker q/.|[l985| ) and for identifying pattern defects (Egolf , 1998 ; Daniels 
\et al. \ |2008| ). Here we applied a two-dimensional Fourier transform to obtain 
the power spectrum within a square, moving window. Local wavelength and 
pattern orientation were identified for each window, and the results averaged 
for all windows that spanned a given location. This straightforward technique 
is suitable for identifying the dominant pattern properties in noisy images with 
irregular patterning. More elegant and rigorously-local techniques, based on the 



Proc. R. Soc. A 6-[29| doi: 10.1098/rspa.00000000 
This journal is © 2011 The Royal Society 



ratios of the spatial derivative of the banding pattern (Egolf, 1998), could not 



be applied because the vegetation bands deviate substantially from a sinusoidal 
pattern in space. The drawbacks of the short-time Fourier transform — namely the 
tendency to truncate long wavelengths due to the finite window size, and poor 
localization of short- wavelength components (Pinnegar & Mansinha, 2004)- 



do not pose significant difficulties in the current application, provided that the 
window size is greater than the local wavelength of the pattern. 

For each window of size L x L, we obtained the 2D fast Fourier transform 
f(k x ,k y ) of the pattern f(x,y). For each window's f(k), we calculated the 

power spectrum S(k) = \f(k)\ 2 . As L increases, the likelihood of identifying 
a single wavelength and orientation decreases. Conversely, as L decreases, the 
^-resolution in Fourier space is reduced. To optimize both the localization of 
measurements and the resolution, we chose L to be at least 5X. The window was 
applied to overlapping regions of the pattern at 20 m intervals. 

The power spectrum measures the power contributed to the pattern by 
each wavevector k. To separate the local wavelength from its orientation, we 



decomposed each k into its magnitude (wavenumber) k = yjk% +kj = 2k /X and 
its orientation 9. Due to symmetry, the pattern orientation can be defined only 
between and K. To identify the dominant k in each window, we binned S(k) 
into annul ar rings of width k = 6n/L.To deconvolve the nat ural \jk scaling of 
the image ( |Burton & Moorhead||1987[|Tolhurst et al\\\992) , we computed the 



total power within each annular ring, S(k). The location peak of this total power 
(rather than the mean power) is used to define the most energetic wavenumber, 
k\. To compensate for the large bin width, we computed the location of the 
weighted average k\ = (kS(k)) / (S(k)) over all rings that formed part of the 
peak and contained > 75% of the peak power. We discarded windows where 
k\ corresponded to the limits of the window function or pixel resolution, and, 
to avoid discarding sites where L truncated A, visually confirmed that these 
windows corresponded to un-patterned areas. Finally, we discarded windows 
where the mean power of the pattern-forming wavenumber was less than that 
of noise. To determine the dominant pattern orientation, we binned S(k) into 
segments of width 7r/8 and computed the average power, S(9). The dominant 
orientation was located via the weighted average 0\ = (OS (6)) / (S(6)) . The most 
energetic values {k\ , 6\ ) for each window were used to generate spatial maps 
A(jc, j) and ®{x,y) representing the local pattern wavelength and orientation. For 
each x,y location, the assigned value of X and comes from the average 2k /k 
and 9 of all windows which include that x,y. 

Although the procedure so far assumes a two-dimensional pattern with 
a single wavelength and orientation within each window, the power spectra 
regularly contained additional peaks. We evaluated the uniqueness of the local 
pattern properties by comparing the dominant peak to its most distant energetic 
peak. Energetic peaks were defined as those containing > 75% of the power 
of the dominant peak. The most distant peak was then defined as the energetic 
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peak with wavenumber k 2 and orientation 02 that maximized \k\ — k 2 \ and 
1 0i — 02 1 - Uniqueness metrics were then defined for both the wavenumber and 
the orientation as: 



Qk = l 



\k\ -k 2 



(2.1) 



Qe = l- 



max(k\ 1 k 2 ) 

I 01 - 02 1 

tt/2 



Qk and Q@ quantify the degree to which the dominant peak is either a unique 
250 energetic peak (£1 = k 2 and 0\ = 02 implies Qk = 1 and Qq = 1) or is one of at 
least 2 energetic peaks lying orthogonal to each other (Qq = 0) or separated by a 
large relative distance in k space. 

Matlab code that performs the Fourier Windowing Method outlined in this 
section is provided as online supplementary material. 



The local Fourier analysis was run twice, using windows with L = 260 m 
and L = 400 m, applying the window to overlapping regions on 20 m intervals, 
and generating smoothly-varying maps of the local pattern properties. The 
two window sizes ensured that the largest wavelengths could be analyzed 

260 (L = 400 m), and exploited the greatest feasible localization given the median 
pattern wavelengths (L = 260 m). The L = 400m method identified large X in 
regions where the L = 260m method identified no pattern; while the L = 260m 
method allowed the pattern properties to be identified in regions with dimensions 
< 400 m. To combine the outcome, we averaged the results of the analyses for 

265 overlapping cells (resulting in a minimum change in the identified and A), and 
retained the uniquely identified and X in cells where either mechanism failed. 
Since the pattern properties in windows located near the edges of the pattern are 
influenced by both patterned and unpatterned regions, we trimmed the resulting 
A (jc,y) and ®{x,y) fields by 130 m to discard the most strongly-affected regions. 

270 The results of the complete analysis are shown in Fig. [2] 

We developed a map of the soil type on the 20 m grid by interpolating the 
SSURGO data. We smoothed the NED elevation data to remove high-resolution 
artifacts ( |Oimoen[ [2000), and computed the topographic gradient (steepness) 
and aspect (orientation) over the same 20 m grid. We made direct comparisons 

275 between X(x,y), ®(x,y) and the soil and topographic features. To cope with the 
large, noisy dataset, we grouped the data into bins of equal size and related the 
median and interquartile range of the local pattern characteristics to the predictor 
variables across the bins. The behavior of these summary statistics and the 
dataset was analyzed using least squares regression. We analyzed the frequency 

280 of occurrence of deviations between the pattern and hillslope orientations A0 and 
the frequency of occurrence of non-unique pattern orientation and wavelength 
Qx and <2© as a function of location (e.g. near topographic minima (drainage 



255 



(c) Construction of the spatial fields used for analysis 
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lines), maxima (ridges), anthropogenic features (roads and trails), and the edges 
of the pattern). 



(d) Comparison to models 

To quantify the degree to which local pattern variations can be explained 
as non-ideal pattern features that can arise from the nonlinear dynamics of the 
pattern forming mechanisms, we performed a control analysis using simulated 
data. Our test region is a 500 x 800 m 2 region of the 0.5 m resolution pattern. 
Analyzed at the 200 x 200 m 2 scale, the pattern in this region contained a 
single high-energy wavenumber and direction peak, indicating that it could be 
reasonably modeled with homogeneous model parameters. We used the pattern 
in this region as the initial condition for a physically-based vegetation band- 
forming model ( [Rietkerk et~aL\ |2002t |Thompson & Katul[ [2009| ). We used 



this model because it (a) represents the sur face runoff feedback me chanisms 
that were observed to occur at this site by |McDonald et al\ ([2008 ); and (b) 



simulates vegetation bands that retain curvature, variability and other non-ideal 
behaviors (compared to more idealized models ( |Lefever & Lejeune} |1997| )) 
and thus has the potential to generate spatially variable pattern properties 
endogenously. We selected model parameters by stepping through reasonable 
parameter combinations and selecting the combination that minimized the 
variation betwee n the observed veg e tation pattern and the model prediction after 
5000 timesteps ( [Thompson et al\ [2008). This method allows us to identify 



a parameter set for the model that approximates the observed patterns as a 
stable steady state solution. Using both the model output and the original binary 
images, we apply the same Fourier windowing method and compute key spatial 
statistics (mean, variance, and autocorrelation length, taken as the lag at which 
the autocorrelation halved) for the resulting X(x,y) and ®(x,y). 



3. Results 

310 (a) Local properties of vegetation patterns 

Using the Fourier windowing method, we find that 44% of the 188 
km 2 area contains patterned vegetation with a clearly identifiable wavelength 
and orientation. The Fourier windowing method fails to detect some areas 
where vegetation patterns are visually identifiable. These include regions where 

315 vegetation bands are confined to fingers of a particular soil type that are < 260 m 
in extent, meaning that the pattern cannot be identified as the dominant Fourier 
mode in any given window. There are also some regions where a pattern can 
be identified by eye, but is so disordered that it falls below the noise threshold. 
Approximately 10% of the study area consists of isolated patches of patterning 

320 which offer little opportunity to explore spatial variations. Instead, we confine 
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our analysis to the 34% 64 km 2 ) of the image that consists of spatially- 
connected vegetation patterns. Within these regions, the pattern has a mean 
wavelength of X = 63 ± 14 m (reported variations are standard deviation unless 
otherwise specified) and an autocorrelation length of approximately 800 m (i.e. 
325 « 12A). The probability distribution function (PDF) of ®(x,y) values is peaked 
in the north-south direction, but spans a full < < n range. Table [T] provides 
summary statistics describing the average pattern characteristics, topography and 
soil properties. 

While seven distinct soil types occur in the study area, two of these, the 

330 Delnorte and Reakor associations, contain 97% of the vegetation patterning 
These soils are distinguished by a relatively low clay and high silt content. The 
Delnorte association is characterized by a very shallow soil depth (23 cm) due 
to the presence of calcium-carbonate based hardpan or petrocalcic horizon. By 
contrast, Reakor association soils are at least 2 m deep. The full range of soil 

335 properties is shown in Table [2j 

We find that three factors broadly determine the local pattern properties: 
the orientation of the slope, the steepness of the slope, and the soil type. The 
orientation of the pattern is almost completely determined by the underlying 
slope orientation. As illustrated in Fig. [3^, approximately 81% of the pattern is 

340 oriented within ±7r/8 radians of the hillslope. Deviations betwee n th e pattern 
and hillslope orientations, denoted A0, are discussed further in § |3|c| Second, 
we observe a significant relationship between the local pattern wavelength X 
and the hillslope gradient. Fig. [3]) illustrates the moderately strong, significant 
decline in the median wavelength X of equally sized data bins with increasing 

345 hillslope gradient (r 2 = 0.63, p = 6 x 10~ 3 ). In addition, we found that the pattern 
wavelength is influenced by the soil type. Fig. ^jp shows PDFs of X for each soil 
type. While Delnorte soils have a unimodal PDF (mode X = 53 m), the Reakor 
association soils have a bimodal distribution. One mode corresponds to that of 
the Deltnorte soils (X = 52 m), but the other modal wavelen gth i s much longer, 

350 with X = 71 m. Further analysis of this effect is provided in §|3|d| 



(b) Comparison to models 

The results are shown in Table [3] We find that while the mean properties 
of the model and observed pattern are similar (a consequence of calibrating 
the model to these means), and X vary 3 — 5 times more in the observed 

355 pattern than in the modeled pattern. The underlying topography is still more 
variable, and presumably causes the variability of the observed pattern. Thus, 
even in the region of the study site with the greatest uniformity in the pattern 
properties, the spatial heterogeneity observed in the real patterns exceeded the 
pattern variability that could be produced by a physical model. These results 

seo justify our attribution of the remaining variability in the pattern properties across 
the site to environmental variation rather than to defects or initial condition 
effects arising from the nonlinear dynamics of the system. 
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(c) Non-uniqueness in pattern attributes 

While the local pattern wavelength and orientation is on average set by the 

365 local hillslope gradient and orientation, we nonetheless observe regions with 
significant deviation from the overall trend. As shown in Fig. [3^, the deviation 
between pattern and hillslope orientation, A©, is > 7r/8 over approximately 20% 
of the pattern. There are even regions in which A0 « 7r/2, where the vegetated 
bands run parallel rather than perpendicular to the local hillslope gradient. Fig. [4] 

370 shows the spatial distribution of the orientation uniqueness metric Qq across 
the patterned region, highlighting regions of non-uniqueness. About half of the 
windows with large A© also contain more than one pattern direction (non- 
uniqueness in ®(x,y)). As illustrated by Fig.|5j non-unique pattern orientations 
are clustered near streams (50% increase in frequency of Qq < 0.75 relative to 

375 the remainder of the pattern), ridges (50% increase in frequency), roads (30% 
increase in frequency), and the pattern edge. We also observe regions with large 
deviations (A© > 7r/8) but nonetheless unique pattern orientations (Qq > 0.75). 
Such regions also occur near ridges and streams more frequently than in the 
remainder of the pattern 30% more frequent in each case). 

380 The insets in Fig. [4] illustrate regions of complexity in (x,y). Fig. [4^ and 
Fig. [4]} show regions of non-unique pattern directions. In Fig. [4^, we show a 
perturbation of the local pattern by a road. As illustrated, the upslope edge 
of roads is globally associated with increased vegetation, while the downslope 
edges are typically bare. Roads thus generate linear features that can create 

385 non-uniqueness in the local pattern orientation. This panel also shows a second 
anthropogenic feature, a storm drainage outlet, which further decreases <2© west 
of the road. Fig. [4J3 illustrates how changes in soil type can lead to low <2©. 
In this example, fingers of the Reakor association soils are interleaved with 
the Delnorte association soils, causing rapid changes in pattern wavelength and 

390 orientation, and consequently multiple energetic values of and X within the 
260 m windows. Fig. [4]: illustrates a region where the pattern orientation changes 
rapidly, turning through approximately K radians within a single 260 m window. 
Such rapid change inevitably leads to low Q® because there are multiple pattern 
orientations located in a single window. There are also, however, regions near 

395 the ridge crest where Q% « 1 and the pattern is locally oriented perpendicular 
to the slope. Fig. |4]1 is representative of the complex vegetation transitions that 
occur near stream channels. Here, too, the pattern lies perpendicular to the local 
hillslope orientation, and rather than curving into the streamline, remains broadly 
aligned with the band orientation away from the stream. 

400 (d) Soil type effects 

While the vegetation patterns on the Delnorte soils have a unimodal 
wavelength distribution, the Reakor soils exhibit a bimodal distribution (see 
Fig. [3]). The pattern associated with the second peak of X on the Reakor soil 
consists of bands of bare soil within a matrix of vegetation cover, inverting the 
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405 distribution in the remainder of the pattern. Visually, this 'anti-stripe' pattern is 
more disordered than the remainder of the pattern seen on the Reakor association 
soils. In the power spectra, these anti-stripe peaks contain less than half the power 
of the short wavelength, coherent peaks. Thus, on the Reakor association soils the 
pattern varies through space from short to long wavelength, ordered to disordered 

410 patterns, and lower to higher biomass. Because increased biomass in drylands 
implies an increased access to water, we examine how X varies as a function 
of the distance to the nearest stream (as a proxy for availability of water). The 
results are shown in Fig. [6] While patterns on the Delnorte soils did not show a 
distance-dependent X (r 2 = 0.08 and p = 0.4), the patterns on Reakor soils show 

415 a trend of increasing X with decreasing distance to the nearest stream (r 2 = 0.96 
and p < 5 x 10" 7 ). 



4. Discussion 

The local Fourier metrics confirm Hypothesis (1) by showing that the wavelength 
and orientation of the vegetation patterns are typically coherent on scales of 

420 600 — 800 m, but can change on scales of « 20 m when the slope orientation 
changes over the same scale. We were unable to reproduce observed variation 
with a model even within regions of relatively uniform patterning, confirming 
Hypothesis (2) and further suggesting that exogenous factors rather than the 
pattern forming mechanisms drive the spatial variability in X and 0. In light 

425 of these findings, and to address Hypothesis (3), our subsequent analysis focuses 
on the relationship between exogenous factors and the local pattern properties. 

Our analysis indicates a hierarchy of controls on the morphology of the 
vegetation patterning: at a global scale the pattern morphology is determined by 
the slope orientation and the hillslope gradient, while the soil type imposes the 

430 template for the pattern forming region. At smaller scales, complexity in the form 
of multiple local pattern orientations and deviations from the global trends are 
associated with ridges, streams, roads and changes in the soil type. Organization 
of the pattern characteristics on small scales arises due to a combination of 
soil type and topographic context: homogeneous, well defined patterns with a 

435 unimodal X on Delnorte soils, regardless of the topographic context. On Reakor 
association soils, the pattern X increases in locations closer to streams, while the 
pattern simultaneously becomes more disordered and vegetation cover increases. 

Several of these relationships support previous observations. The average 
hillslope gradient of 0.7% at the Fort Stockton site conforms to observations and 

440 predictions of a minimum hillslope gradient being needed for band formation 
( [Ursino & Contarini]|2006[ ), and exceed s the minimum slope gra dient observed in 
other settings (e.g. 0.25% in the Sudan (Debla uwe et fl/.[|201 In addition, the 
inverse relationship between pattern wavelength and the hillslope gradient (see 



Fig. [3]) has been identified elsewhere ([Valentin et al\ 


1999;d'Herbes etaL, 2001 


Eddy et al 1999 ; Deblauwe et al 


201 1|). However, the one model (JKlausmeier 
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1999 ) that has been used to study wavelength-gradien t correla t ions predicts the 
opposite trend: steepness causes larger wavelengths (She rratt[ |2005[ |Ursino & 



Contarini, 2006 ). We explored this trend with a different model ( fRietkerk et ~al. 



2002]) and find that it also produces the opposite trend to that observed in natural 
systems. While both models predict that multiple pattern wavelengths are stabl e 
for a given hillslope gradient ( Sherrat t~& Lord||2007t |T"hompson & Katul, 2009 ), 



the selection of particular longer or shorter wavelengths within that range is 
clearly at odds with observations. 

The close correspondence between hillslope and pattern orientation at the 
Fort Stockton site is a near-universal feature of vegetation patterning. The 
observed deviations generally arise due to local effects that disrupt the pattern, 
rapidly alter the direction of water flow, or change the strength and/or length- 
scales of the pattern forming feedbacks. 

Pattern disruption is exemplified by the effects of roads (see Fig. [4j 
Roads disconnect upslope and downslope vegetation bands by preventing the 
redistribution of water. Rapid alteration of the direction of water flow occurs 
along ridges and stream channels. The pattern near these locations is less likely to 
have a unique orientation or wavelength than in the mid-slope areas, and is more 
likely to have a large A®. The prominent ridgeline shown in Fig. |4j: is locally 
surrounded by a region with A© « k/2 and Qq > 0.75. This location provides an 
unambiguous example of a change in pattern ori entation near no-flow boundary 
conditions, as predicted by McGra thgf q/.| ( |2012| ). Other ridgelines are associated 



with large A©, but the low Q% in these locations makes the interpretation of the 
observed A© ambiguous. Methods that identify a truly local metric of pattern 
properties, instead of the quasi-local metric used here could both help to resolve 
these ambiguities, and to extend the analysis into patterned areas less than 260 m 
wide. 

Stream channel locations are also associated with rapid, and sometimes 
discontinuous changes in pattern orientation. Like the ridge shown in Fig. [4b, 
the stream channel shown in Fig.[4]l is surrounded by a region where A© « k/2 
and Qq > 0.75, i.e. a unique pattern orientation that is perpendicular to the 
one expected from the slope orientation. There has been little exploration of 
the interaction of vegetation pattern with stream channels, presumably because 
the current paradigm of vegetation pattern models offer little reason to think 
that such interactions would be important. Stream channels occur in locations 
where surface runoff rapidly flows away and therefore cannot affect vegetation 
upslope from the channel. However, we find evidence that the distance from a 
stream channel alters the pattern properties on the Reakor Association soil type. 
On these soils the pattern wavelength, vegetation cover and disorder increase 
near the streams. These observations suggest that an additional mechanism could 
affect the patterning at this study site: for instance, the stream-channel boundary 
condition propagating back up into the hillslope. 

We hypothesize that the ultimate cause of the changes in pattern morphology 
on the Reakor soils lies in the contrast in the soil depth between the Delnorte and 
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Reakor associations: from 23 cm to over 2 m. We are not the first to propose an 
association between soil depth and the vegetation pattern structure. Depth to a 
silcrete hardpan is associated with a transition between sharp patterns (shallow 



h 
F 



ardpan) and diffuse patterns (deep or no hardpan) in Australia (Mabb utt & 
Fanning[[1987[|Tongway & Ludwig [T990). Strikingly cohere nt veg etation bands 



in Niger occur above a shallow ironstone har dpan (| White) |1970|). Th e broad, 



diffuse bandin g near Fort Stockton st udied by McDonal d et al\ ( |2008| ) occurs 
on deep clays. |McDonald et <z/.| ( |2008] ) cited unpublished research claiming that 



vegetated bands were associated with a local increase in the depth to the hardpan, 
similar to previous observations of increased vegetation density on deep soils in 
Australia ( [Mott & McCombj [T974] ). 



Field evidence is needed to determine the exact mechanisms by which soil 
depth and proximity to streams could lead to the changes we observed in pattern 
morphology. Three scenarios illustrating potential mechanisms are shown in 
Fig. [7] and could form the basis for future field studies. First, shallow soils 
could promo te lateral root extension by plants , exaggerating the effects of root 
competition ( [Gilad et at] |2004[ |Yizhaq et at] |2005| ) (see Fig. [TJib). Studies of 



Chihuahuan desert species confirm that there is intense root competition in the 



^one abo ve the petrocalcic horizon, where root growth is concentrated fGibbens 
|& Lenzl|200TT ). 



In the second scenario, we recognize that shallow surface soils have limited 
water storage and saturate readily. For example, the Delnorte Association has as 
little as 2 cm of water storage capacity in the soils above the impeding layer. 
Unlike dry soils, which only generate runoff during very heavy rains, saturated 
soils shed all rainfall as runoff. If soils did not saturate near plants this runoff 
water could be trapped and infiltrate in these locations (see F ig. [7b) . Plant roots 
penetrate and thus may break up hardpans (Gibbe ns~& Lenz[|2 001 ). Root water 



uptake is also a driver of hardpan formation (Duni way fl/.[|2007| ), and hardpans 
might thus form at greater depth beneath deep-rooted (woody) vegetation. Either 
mechanism could prevent the surface soil from saturating near the vegetated 
bands and allow it to store runoff. 

The third scenario also relates to the potential for saturation to occur above 
the hardpan, allowing subsurface saturated flow to occur (see Fig. [7]:). The 
changes in pattern properties with distance to the stream suggests that water 
availability increases downslope: this requires subsurface storage, if not flow. A 
relationship between the pattern formation and subsurface flow could explain the 
local deviations of pattern orientation from slope orientation near the streamlines 
since the surface orientation might not correspond to the local water table flow 
direction. 

All three scenarios would tend to strengthen the pattern forming feedbacks 
on the shallow soils. They could be investigated using subsurface soil moisture 
sensors to observe shallow soil saturation; water isotope tracers to determine 
the water sources used by plants, observations of calcium ion concentrations in 
runoff which would provide an indicator of water contact with the petrocalcic 
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horizon and root excavations to compare morphologies in sites with different 
depths to the impeding layer. These studies could be valuable to provide more 
information about the hydrological role of petrocalcic horizons, which have 
received relatively little attention given their ubiquity in desert environments 
( [Duniway efaT^Ml) . 



5. Conclusion 

Large scale analyses of variation in pattern morphology has provided broad 
confirmation of many theoretical predictions about vegetation patterning and its 
variation along climatic gradients (Debla uwe et fl/.||2011| ), which are consistent 
with our observations at local scales. Our analysis shows that while pattern 
properties mostly vary on scales of 600 — 800 m, they can change much more 
rapidly around boundaries in topography or soil characteristics. By analyzing 
the pattern on these fine scales, we identified deviations between hillslope and 
pattern orientation, soil-controlled changes in pattern wavelength, coherence and 
vegetation cover, and at least one like l y exa mple of the boundary condition 
effects predicted by |McGrath et al.\ ( |2012| ). Two observations, echoed at 



multiple other sites, are not well-explained by current theory: the observation of 
increasing pattern wavelength on steepness, and the soil controls on vegetation 
pattern length-scale and coherence. 

The large timescale separation between individual storm events and the 
timescales on which desert vegetation distributions change means that an 
ongoing dialogue between empirical and theoretical studies is critical for 
understanding the dynamics of these ecosystems. Local information about 
pattern qualities, when combined with high resolution information about the 
pattern substrate, is evidently a useful additional tool for analysis. 

Despite the advances in understanding vegetation patterns in the past 
ten years, theoretical models still require the use of effective parameters to 
describe feedbacks, soil-plant-water interactions and the resulting landscape 
fluxes. Linking theory and observation to make quantitative predictions, 
therefore, remains an outstanding challenge. Addressing this challenge requires 
improved observations of within-storm hydrologic processes and plant water 
use: observations that ecohydrologists are increasingly equipped to make 
due to developments in distributed sensing systems and tracer technologies. 
Computational tools for assessing three-dimensional soil moisture dynamics, 
land-atmosphere interactions and vegetation spread are also improving. By 
coupling these tools with detailed field measurements, there is potential to 
develop a detailed theoretical framework that can address the consequences 
of changing soil depth, root orientation, runoff generation mechanisms and 
subsurface flow processes on the overall dynamics and resilience of vegetation 
patterns. 
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7. Tables 



Table 1 . Pattern characteristics of natural patterning across the full pattern extent. 





Mean 


Std. Dev. 


Corr. Length (m) 


Pattern wavelength, A (m) 


63 


13.8 


800 


Mean pattern orientation, (rads) 


1.4 


0.76 


600 


Hillslope orientation (rad) 


1.2 


0.75 


1400 


Hillslope gradient 


0.007 


0.002 


1900 



Table 2. Soil properties for the Delnorte and Reakor Association Soils. Note that the soil textural 
percentages do not sum to 100 because gravel and organic content have not been reported. 





Reakor 


Delnorte 


Soil depth (cm) 


200 


23 


Hydraulic conductivity (K sat , jum s -1 ) 


2.1 


59.2 


Clay content (%) 


31.5 


8.7 


Sand content (%) 


6.8 


42.8 


Silt content (%) 


61.7 


21.1 



Table 3. Pattern characteristics comparing analysis of natural and simulated data for a a 500 x 
800 m 2 region. The simulated data comes from a modified Rietkirk model calibrated to match the 
mean pattern properties of the natural data. 





Mean 


Std. Dev. 


Corr. Length (m) 


Natural pattern wavelength A (m) 


53 


4.1 


350 


Model pattern wavelength A. (m) 


54 


1.7 


360 


Natural pattern orientation (rad) 


1.7 


0.39 


420 


Model pattern orientation (rad) 


1.6 


0.07 


390 


Hillslope orientation (rad) 


1.7 


0.5 


120 



8. Figures 
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Figure 2. Maps of (a) local wavelength k(x,y), with color scale indicating wavelength in meters 
and (b) local orientation ©(x,y), with color scale indicating orientation in radians. Areas in white 
do not contain patterns as recognized by the Fourier windowing method. 




Figure 3. (a) Frequency distribution of A®, the deviation between the local pattern orientation 
and the local topographic orientation. The slope orientations deviate by less than 7r/8 for 80% of 
the image, (b) Comparison of the local pattern wavelength A to the local topographic slope. Data 
were collected into 10 bins of equal size 1 1000 data points per bin). The box and whisker plots 
indicate the median (central lines); the 25th and 75th quartiles (box limits) within the bin. The 
whiskers extend two inter-quartile ranges from the median. The thick red line indicates a linear 
fit to the medians, (c) PDF of the pattern wavelengths associated with the two major soil classes, 
indicating the bimodality in A on Reakor association soils. 
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Figure 4. Spatial distribution of the uniqueness metric <2® f° r the pattern orientation, and 
magnified cases illustrating different examples of non-uniqueness in the pattern orientation. Panel 
(a) indicates non-uniqueness in orientation associated with a road, which adds north-south oriented 
structure to the local NE-SW oriented banding due to upslope ponding of runoff. Panel (b) 
indicates a region of changing pattern orientation and wavelength associated with the interleaving 
of two contrasting soil types. Panel (c) illustrates the short length-scales on which the pattern can 
change orientation as it bends around a sharp ridge (location indicated in blue). A proportion of 
the region where A® « tf/4, as indicated by the different director arrows, has Qq > 0.75. Panel (d) 
indicates the deviation between pattern and slope orientation around a stream channel (estimated 
location shown in blue); again there are regions where A@ « n/2 and Q® > 0.75. 
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Figure 5. Relative increase in frequency of (a,b) non-unique pattern orientations and (c) 
wavelengths, for windows within 250 m of local elevation minima (streams), maxima (ridges), 
roads, and pattern boundaries. All frequencies are measured relative to the frequency in the 
remainder of the pattern. Non-unique orientations, for example, occurred 50-60% more often near 
streamlines and ridges than in the bulk of the pattern. 
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Figure 6. Distribution of the soils within the pattern forming areas with respect to landscape 
position, referenced to the nearest stream channel or local elevation minimum for (a) Reakor 
and (b) Delnorte Association soils. Although both soil types occur across the range of 
landscape positions, the Reakor association soils are more strongly associated with riparian areas. 
Distribution of pattern wavelengths on the Reakor Association (c) and Delnorte Association (d) 
soils with respect to the nearest stream, again based on equally-sized bins of the dataset. There is a 
strong (r 2 = 0.96) and significant (p < 5 x 10 -7 ) decline in the wavelength on Reakor soils when 
moving from the riparian areas to the uplands. There is no correlation between landscape position 
and pattern properties on the Delnorte association (r 2 = 0.08 and p = 0.4). 
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Figure 7. Schematic illustration of possible effects of a shallow impeding layer on pattern 
formation mechanisms, (a) In the absence of an impeding layer, plants develop deep root systems 
to exploit water from across the whole soil profile (b) A shallow impeding layer restricts the 
vertical growth of roots, confining them to the surface soils, and promoting root competition 
laterally, (c) The shallow impeding layer also creates the potential for complex subsurface 
hydrology. Saturation of soils above the impeding layer may promote saturation excess runoff. 
Gradients - due to changes in the depth of the impeding layer, or in water potential across the 
impeding layer - can also induce saturated porous media flow. Provided that surface soils near 
plants can freely drain, the altered hydrology may exaggerate the positive feedback between plants 
and water. 
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